/*------------------------------------------------------------------------
 *  Copyright (C) 2011  Luis M. de la Cruz Salas
 *
 *  This file is part of TUNA
 *
 *  TUNA is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  TUNA is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ------------------------------------------------------------------------*/

#ifndef _NONUNIFORM_H_
#define _NONUNIFORM_H_

#include "Meshes/StructuredMesh.hpp"
#include <iostream>
#include <fstream>
#include <string>
#include <vector>

namespace Tuna {

  template <typename T, int Dim> class NonUniform;

/*!
 *  Specialization for uniform meshes in 1D.  
 */
  template <typename T>
  struct Typeinfo <NonUniform<T, 1> > {
    typedef T prec_t;
    enum { Dim = 1 };
  };

/*!
 *  Specialization for uniform meshes in 2D.  
 */
  template <typename T>
  struct Typeinfo <NonUniform<T, 2> > {
    typedef T prec_t;
    enum { Dim = 2 };
  };

/*!
 *  Specialization for uniform meshes in 3D.  
 */
  template <typename T>
  struct Typeinfo <NonUniform<T, 3> > {
    typedef T prec_t;
    enum { Dim = 3 };
  };

/*! 
 ***************************************************************************
 *  Adaptor for NON-uniform structured meshes on rectangular domains.
 *  \author  Luis M. de la Cruz [ Sun Jan 20 14:16:36 GMT 2010 ]
 ***************************************************************************
 */
template<class Tprec, int Dim>
class NonUniform : public StructuredMesh<NonUniform<Tprec, Dim> >  {
  
// Variables from StructuredMesh classes
  using StructuredMesh< NonUniform< Tprec, Dim > >::xyz; // this comes from Mesh
  using StructuredMesh< NonUniform< Tprec, Dim > >::dxyz;
  using StructuredMesh< NonUniform< Tprec, Dim > >::dface;

 public:

  NonUniform() : StructuredMesh<NonUniform<Tprec, Dim> >() { }
  ~NonUniform() {}

  bool print() {
    std::cout << " NonUniform < " << Dim << "D > >"
	      << "\n +-----+";
  }

  template<typename Tfunct>
  inline bool calcPoints(Axis_t, Tfunct, Tprec, Tprec, int);
  inline bool calcDeltas(Axis_t);
protected:    

private:

};
  
/*!
 *  This function calculates the distribution of the points in the
 *  axis defined by the \c a parameter using the function \c f. The
 *  function is applied between \c x_1 and \c x_2 .
 *  \param a axis
 *  \param f functor that defines the distribution
 *  \param x_1 begining of the segment
 *  \param x_2 finish of the segment
 *  \param num_nodes total number of number used in the segment
 */
template<typename Tprec, int Dim>
template<typename Tfunct>
inline bool
NonUniform<Tprec, Dim>::calcPoints(Axis_t a, Tfunct f, 
				   Tprec x_1, Tprec x_2, 
				   int num_nodes) 
{
  Tprec delta = 1.0 / (num_nodes - 1); // in the computational space
  Tprec fact = x_2 - x_1;              // in the physical space

// this->... refers to variables defined in the Mesh class. 
  xyz[a].pop_back(); // first point  
  for(int i = 0; i < num_nodes; ++i) 
    xyz[a].push_back ( x_1 + fact * f(i*delta) ); // apply the function

  return 0;
}

template<typename Tprec, int Dim>
inline bool
NonUniform<Tprec, Dim>::calcDeltas(Axis_t a)
{
  //       x(i-1)   x(i)    x(i+1)          <--- nodes numeration 
  //   |     |   i   |   i+1   |           |    
  //   o--*--|---*---|----*----|-----*-----o    * - internal volumes
  //   |     |<----->|<------->|           |    o - boundary volumes
  //             /        \                 \    
  //     dxyz[a][i]     dxyz[a][i+1]         NOTE: the size of boundary                 
  //                                                 volumes is zero
  int N = xyz[a].size();      // Number of nodes
  dxyz[a].push_back( 0.0 );   // First delta is zero
  for(int i = 1; i < N; ++i)
    dxyz[a].push_back( xyz[a][i] - xyz[a][i-1] );
  dxyz[a].push_back( 0.0 );    // Last delta is zero


  //       x(i-1)   x(i)    x(i+1)               
  //   |     |   i   |   i+1   |           |    
  //   o--*--|---*---|----*----|-----*-----o   
  //      |      |        |
  //      |<---->|<------>|   
  //          /       \                 
  //   dface[i-1]    dface[i]  
  // 
  // First and last deltas should be multiplied by 0.5, but this factor
  // is (should be) taken into account in boundary implementations.
  //
  dface[a].push_back( (xyz[a][1] - xyz[a][0]) * 1.0 );     // First delta
  for(int i = 1; i < N - 1; ++i)              
    dface[a].push_back( ( xyz[a][i+1] - xyz[a][i-1] ) * 0.5 );
  dface[a].push_back( (xyz[a][N-1] - xyz[a][N-2]) * 1.0 ); // Last delta

  return 0;
}

} //  Tuna namespace
 
#endif // _NONUNIFORM_H_








